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Abstract. We present results of numerical computations of quasiequilibrium 
1 sequences of binary neutron stars with zero vorticity, in the general relativistic 

' framework. The Einstein equations are solved under the assumption of a conformally 

^sO ' flat spatial 3-metric (Wilson-Mathews approximation) . The evolution of the central 

density of each star is monitored as the orbit shrinks in response to gravitational 
wave emission. For a compactification ratio M/R = 0.14, the central density 
remains rather constant (with a slight increase, below 0.1%) before decreasing. For 
^ ' a higher compactification ratio M/R = 0.17 (i.e. stars closer to the maximum 

' mass configuration), a very small density increase (at most 0.3%) is observed before 

' the decrease. This effect remains within the error induced by the conformally flat 

^ approximation. It can be thus concluded that no substantial compression of the stars 

, is found, which would have indicated a tendency to individually collapse to black hole 

— prior to merger. Moreover, no turning point has been found in the binding energy 

or angular momentum along evolutionary sequences, which may indicate that these 
systems do not have any innermost stable circular orbit (ISCO). 
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Introduction 



Inspiraling neutron star binaries are expected to be among the strongest sources 
5^ ■ of gravitational radiation that could be detected by the interferometric detectors 

\ currently under construction (GEO600, LIGO, TAMA and Virgo). These binary 

systems are therefore subject to numerous theoretical studies. Among them are (i) 
Post-Newtonian (PN) analytical treatments (e.g. Q, (D, Q) and (ii) fully relativistic 
hydrodynamical treatments, pioneered by the works of Oohara and Nakamura (see 
e.g. and Wilson et al. (5[ ^. The most recent numerical calculations, those 
of Baumgarte et al. 0, || and Marronetti et al. rely on the approximations 
of (i) a quasiequilibrium state and (ii) of synchronized binaries. Whereas the first 
approximation is well justified before the innermost stable orbit, the second one does 
not correspond to physical situations, since it has been shown that the gravitational- 
radiation driven evolution is too rapid for the viscous forces to synchronize the spin 
of each neutron star with the orbit |l^ as they do for ordinary stellar binaries. 
Rather, the viscosity is negligible and the fluid velocity circulation (with respect to 
some inertial frame) is conserved in these systems. Provided that the initial spins are 
not in the millisecond regime, this means that close binary configurations are better 
approximated by zero vorticity (i.e. irrotationat) states than by synchronized states. 
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Dynamical calculations by Wilson et al. ^, || indicate that the neutron stars may 
individually collapse into a black hole prior to merger. This unexpected result has 
been called into question by a number of authors (see Ref. for a summary of all 
the criticisms and their answers). Recently Flanagan p3[ has found an error in the 
analytical formulation used by Wilson et al. |^. This error may be responsible of 
the observed radial instability. As argued by Mathews et al. jl^, one way to settle 
this crucial point is to perform computations of relativistic irrotational configurations. 
We have performed recently such computations . They show no compression of the 
stars, although the central density decreases much less than in the corotating case. In 
the present report, we give more details about the results presented in Ref. ||l4| and 
extend them to the compactification ratio M/R = 0.17 (the results of Ref. [[l4| have 
been obtained for a compactification ratio M/R = 0.14). 



2. Analytical formulation of the problem 



2.1. Basic assumptions 

We have proposed a relativistic formulation for quasiequilibrium irrotational binaries 
[ psf as a generalization of the Newtonian formulation presented in Ref. The 
method was based on one aspect of irrotational motion, namely the counter-rotation 
(as measured in the co-orbiting frame) of the fluid with respect to the orbital motion 
(see also Ref. [0). Since then, Teukolsky Q and Shibata [|l^ gave two formulations 
based on the definition of irrotationality, which implies that the specific enthalpy 
times the fiuid 4- velocity is the gradient of some scalar field {potential flow) . The 
three formulations are equivalent; however the one given by Teukolsky and by Shibata 
greatly simplifies the problem. Consequently we used it in the present work. 

The irrotational hypothesis amounts to say that the co-momentum density is the 
gradient of a potential 

/iu = V*, (1) 

where h and u are respectively the fluid speciflc enthalpy and fiuid 4-velocity. 

Beside the physical assumption (|l|), two simplifying approximations are 
introduced: 

(i) The spacetime is supposed to have a helicoidal symmetry JTst , which means that 
the orbits are exactly circular and that the gravitational radiation content of 
spacetime is neglected. 

(ii) The spatial 3-metric is assumed to be conformally flat (Wilson-Mathews 
approximation ||2l] , |^), so that the full spacetime metric reads 

ds^ = -{N^ - BiB')dt^ - 2B^dt dx' + A^fijdx' dx^ , (2) 



where is the flat space metric. 

The Killing vector corresponding to hypothesis (i) is denoted by 1 (cf. Fig. |l|). 

Approximation (i) is physically well motivated, at least up to the innermost 
stable orbit. Regarding the second approximation, it is to be noticed that (i) the 
1-Post Newtonian (PN) approximation to Einstein equations fits this form, (ii) it is 
exact for arbitrary relativistic spherical configurations and (iii) it is very accurate 
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Figure 1. Spacetime foliation Et, helicoidal Killing vector I" and its trajectories 
= const, which are the worldlines of the co-orbiting observer (4-velocity: v°'). 
Also shown are the rotating-coordinate shift vector and the unit future-directed 
vector n", normal to the spacelike hypersurface St (Figure from Ref. [ta]). 



for axisymmctric rotating neutron stars p2| . An interesting discussion about some 
justifications of the Wilson-Mathews approximation may be found in A stronger 
justification may be obtained by considering the 2.5-PN metric obtained by Blanchet 



et al. 1 23 for point mass binaries. Using Eq. (7.6) of Ref. |23|, the deviation from a 
conformally flat 3-metric (which occurs at the 2-PN order) can be computed at the 
location of one point mass (i.e. where it is expected to be maximal), the 3-metric hij 
being written as 

h., = A' f., + h^^ + h^f^^ . (3) 

The result is shown in Table |l| for two stars of IAMq each. It appears that at a 
separation as close as 30 km, where the two stars certainly almost touch each other, 
the relative deviation from a conformally flat 3-metric is below 2%. 



2.2. Equations to be solved 

We refer to Ref. jwj for the presentation of the partial differential equations (PDEs) 
which result from the assumptions presented above. In the present report, let us 
simply mention a point which seems to have been missed by various authors: the 
existence of the flrst integral of motion 

hi ■ u = const. (4) 

does not result solely from the existence of the helicoidal Killing vector 1. Indeed, 
Eq. (^ is not merely the relativistic generalization of the Bernoulli theorem which 
states that ft, 1 • u is constant along each single streamline and which results directly 
from the existence of a Killing vector without any hypothesis on the flow. In order for 
the constant to be uniform over the streamlines (i.e. to be a constant over spacetime), 
as in Eq. (^), some additional property of the flow must be required. One well known 
possibility is rigidity (i.e. u colinear to 1) p4| . The alternative property with which 
we are concerned here is irrotationality [Eq. (Ilj)]. This was first pointed out by Carter 
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Table 1. Deviation from a conformally flat 3-metric at 2-PN at 2.5-PN order for 
point mass binaries of M = 1.4 Mq each, d is the coordinate separation between the 
two stars, v is the coordinate orbital velocity of one star, t<J2PN is the orbital angular 
velocity, the other notation s ar e defined by Eq. (^. The metric at the 2.5PN level is 
taken from Blanchet et al. ]23l and is computed at the location of the point masses. 
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Figure 2. Logarithm of the relative global error of the numerical solution with 
respect to the number of degrees of freedom in 9 for a Roche ellipsoid for an equal 
mass binary system and /{irGp) = 0.1147 (the numbers of degrees of freedom in 
the other directions are Nr = 2Ng — 1 and = Ng — I). Also shown is the error in 
the verification of the virial theorem (Figure from Ref. |26| ) . 



3. Numerical procedure 

3.1. Description 

The numerical procedure to solve the PDE system is based on the multi-domain 
spectral method presented in Ref. We simply recall here some basic features of 
the method: 

• Spherical- type coordinates {S,,9',ip') centered on each star are used: this ensures 
a much better description of the stars than with Cartesian coordinates. 
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• These spherical-type coordinates are surface-fitted coordinates: i.e. the surface 
of each star Hes at a constant value of the coordinate ^ thanks to a mapping 
{£,,9',ip') {r,9,(p) (see for details about this mapping). This ensures that 
the spectral method applied in each domain is free from any Gibbs phenomenon. 

• The outermost domain extends up to spatial infinity, thanks to the mapping 
1/r = (1 — ^)/(2i?o)- This enables to put exact boundary conditions on the 
elliptic equations for the metric coefficients: spatial infinity is the only location 
where the metric is known in advance (Minkowski metric). 

• Thanks to the use of a spectral method in each domain, the numerical error 
is evanescent, i.e. it decreases exponentially with the number of coefficients (or 
equivalently grid points) used in the spectral expansions, as shown in Fig. ^. 

The PDE system to be solved being non-linear, we use an iterative procedure. This 
procedure is sketched in Fig. ^. The iteration is stopped when the relative difference 
in the enthalpy field between two successive steps goes below a certain threshold, 
typically 10"'^ (cf. Fig. 



The Equalioii Of Stale 
The central enthalpies HI and H2 
The separation between the stars 



Compute 

The orbital frequency 
The location of the rotation a.K; 



Solve ^ Solve 

m-like field equations for N Poisson-like field equations for B' 



^ Solve ^ Solve 

The Poisson-like equations for I'l The Poisson-like equations for '}'2 



Compute 
The enthalpy distribution 2 



Compute 

The coordinnate mappings 



Figure 3. Schematic representation of the iterative procedure used in the numerical 
code. 



The numerical code is written in the lorene language [Q, which is a C++ based 
language for numerical relativity. A typical run make use of A^^ = 33, Nq = 21, 
and = 20 coefficients {= number of collocation points, which may be seen as 
number of grid points) in each of the domains on the multi-domain spectral method. 
8 domains are used : 3 for each star and 2 domains centered on the intersection between 
the rotation axis and the orbital plane. The corresponding memory requirement is 
155 MB. A computation involves ~ 250 steps (cf. Fig. which takes 9 h 30 min 
on one CPU of a SGI Origin200 computer (MIPS RIOOOO processor at 180 MHz). 
Note that due to the rather small memory requirement, runs can be performed in 
parallel on a multi-processor platform. This especially useful to compute sequences of 
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Convergence curve 

Resu1208a/E7: irrotational configuration IVIg=1.625 IVI, 

10"' , : 




50 100 150 200 250 300 

step 

Figure 4. Convergence (measured by the relative difference in the enthalpy field 
between two successive steps) of the iterative procedure for one of the irrotational 
models considered in Ref. [M. The bump around the 90th step corresponds to the 
switch on of the procedure of convergence toward a given baryon mass. 



configurations. 



3.2. Tests passed by the code 

In the Newtonian and incompressible Hmit, the analytical solution constituted by a 
Roche ellipsoid is recovered with a relative accuracy of ~ 10~^, as shown in Fig. |^. For 
compressible and irrotational Newtonian binaries, no analytical solution is available, 
but the virial theorem can be used to get an estimation of the numerical error: we 
found that the virial theorem is satisfied with a relative accuracy of 10~^. A detailed 
comparison with the irrotational Newtonian configurations recently computed by Uryu 
& Eriguchi (29[ |o| will be presented elsewhere. Regarding the relativistic case, we have 
checked our procedure of resolution of the gravitational field equations by comparison 
with the results of Baumgarte et al. j|] which deal with corotating binaries [our code 
can compute corotating configurations by setting to zero the velocity field of the fluid 
with respect to the co-orbiting observer] . We have performed the comparison with the 
configuration za = 0.20 in Table V of Ref. |^. We used the same equation of state 
(EOS) (polytrope with 7 = 2), same value of the separation rc and same value of 
the maximum density parameter q"^^^. We found a relative discrepancy of 1.1% on f2, 
1.4% on Mo, 1.1% on M, 2.3% on J, 0.8% on za, 0.4% on and 0.07% on rg (using 
the notations of Ref. Q). 



4. Numerical results 
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4-.1. Equation of state and compactification ratio 

As a model for nuclear matter, we consider a polytropic equation of state (EOS) with 
an adiabatic index 7 = 2: 

p = K(mBn)''' , e — rriBn + p/ ["f ~ I) , (5) 

where p, n, e are respectively the fluid pressure, baryon density and proper energy 
density, and me — 1-66 x 10~^^ kg, k = 1.8 x 10~^ Jm'^kg^^. This EOS is the same 
as that used by Mathews, Marronetti and Wilson (Sect. IV A of Ref Q). 

Spherical models in isolation 

Polytrope y=2 k=0.0332 p„,„"'c'^ 

2.2 
2 
1.8 
— 1.6 

CO 

w 

CO 

s 

1.4 
1.2 
1 

2 4 6 8 10 12 14 16 18 
Central energy density e„ [ p^i.^c'' ] 

Figure 5 . Mass as a function of the central energy density for static isolated neutron 
stars constructed with the EOS The heavy dots are configuifltioiis considered 
by our group |14| and Marronetti, Mathews & Wilson group |L3], [B1| (see text) 
(Pnuc := 1.66 X 10" kgm-3). 

The mass - central density curve of static configurations in isolation constructed 
upon this EOS is represented in Fig. ^ The three points on this curve corresponds 
to three configurations studied by our group and that of Marronetti, Mathews and 
Wilson: 

• The configuration of baryon mass A/b — 1.625 Mq and compactification ratio 
M/R — 0.14 is that considered in the dynamical study of Mathews, Marronetti 
and Wilson ||l^ and in the quasiequilibrium study of our group (Ref. [Q and 
this paper). 

• The configuration of baryon mass Mb = 1.85 Af© and compactification ratio 
M/R = 0.17 is studied in the present paper. 

• The configuration of baryon mass A^b = 1.95 Af0 and compactification ratio 
Af/i? = 0.19 has been studied recently by Marronetti, Mathews and Wilson |3l[||t| 
by means of a new code for quasiequilibrium irrotational configurations. 

t Marronetti et al. ^ij use a different value for the EOS constant k: their baryon mass Mb = 
1.55 M0 must be rescaled to our value of k in order to get Mb = 1.95 M©. Apart from this scaling, 
this is the same configuration, i.e. it has the same compactification ratio M/R = 0.19 and its relative 
distance with respect to the maximum mass configuration, as shown in Fig. ral is the same. 
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^.2. Irrotational sequence with M/R — 0.14 

In this section, we give some details about the irrotational sequence Mb — 1.625 M0 
presented in Ref. |lj]. This sequence starts at the coordinate separation d — HQ km 
(orbital frequency f — 82 Hz), where the two stars are almost spherical, and ends at 
d ~ Al km (/ = 332 Hz), where a cusp appears on the surface of the stars, which 
means that the stars start to break. The shape of the surface at this last point is 
shown in Fig. |^. 



Figure 6. Shape of irrotational binary neutron stars of baryon mass Mb = 
1.625 Mq, when the coordinate separation between their centers (density maxima) is 
41 km. Only one half of each star is represented (the part which is above the orbital 
plane). The drawing is that of the numerical grid, which coincides with the surface 
of the star, thanks to the use of surface-fitted spherical coordinates. 

The velocity field with respect to the co-orbiting observer, as defined by Eq. (52) 
of Ref. Il^ , is shown in Fig. 0. Note that this field is tangent to the surface of the 
star, as it must be. 

The lapse function TV (cf. Eq. ||) is represented in Fig. |8[ The coordinate system 
(x, y, z) is centered on the intersection between the rotation axis and the orbital plane. 
The X axis joins the two stellar centers, and the orbital is the z = plane. The value 
of N at the center of each star is = 0.64. 

The conformal factor of the 3-metric [cf. Eq. (||)] is represented in Fig. ||. Its 
value at the center of each star is — 2.20. 

The shift vector of nonrotating coordinates, N, (defined by Eq. (9) of Ref. ||lj]) is 
shown in Fig. |l^. Its maximum value is 0.10 c. 

The Kxy component of the extrinsic curvature tensor of the hypersurfaces t = const 
is shown in Fig. |l^. We chose to represent the Kxy component because it is the only 
component of Kij for which none of the sections in the three planes x = Q, y = Q and 
2 = vanishes. 

The variation of the central density along the M-q = 1.625 Mq sequence is shown 
in Fig. £2l We have also computed a corotating sequence for comparison (dashed line 
in Fig. 03). In the corotating case, the central density decreases quite substantially as 
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Velocity / co — orbiting observer 

I , , , , I , , , , I 




-10 10 

X [l^m] 

Figure 7. Velocity field with respect to the co-orbiting observer, for the configuration 
shown in Fig. a. The plane of the figure is the orbital plane. The heavy line denotes 
the surface of the star. The companion is located at x = +41 km. 



the two stars approach each other. This is in agreement with the resuhs of Baumgarte 
et ah 0, ||. In the irrotational case (sohd hne in Fig. |l^), the central density remains 
rather constant (with a shght increase, below 0.1%) before decreasing. We can thus 
conclude that no tendency to individual gravitational collapse is found in this case. 
This contrasts with results of dynamical calculations by Mathews et al. which show 
a central density increase of 14% for the same compactification ratio M/R = 0.14. 



4-.3. Irrotational sequence with M/R — 0.17 

In order to investigate how the above result depends on the compactness of the stars, 
we have computed an irrotational sequence with a baryon mass Afe = 1.85 Mq, which 
corresponds to a compactification ratio M/R — 0.17 for stars at infinite separation 
(second heavy dot in Fig. ||). The result is compared with that of M/R = 0.14 in 
Fig. O. A very small density increase (at most 0.3%) is observed before the decrease. 



Note that this density increase remains within the expected error 2%, cf. Sect. 2.1) 
induced by the conformally flat approximation for the 3-metric, so that it cannot be 
asserted that this effect would remain in a complete calculation. 

Marronetti, Mathews and Wilson [ pTf have recently computed quasiequilibrium 
irrotational configurations by means of a new code. They use a higher compactification 
ratio, M/R = 0.19 (third heavy dot in Fig. ||). They found a central density increase 
as the orbit shrinks much pronounced than that we found for the compactification 
ratio M/R = 0.17: 3.5% against 0.3%. We will present irrotational sequences with 
the compactification ratio M/R = 0.19 and compare with the results by Marronetti 
et al. ISlll in a future article. 
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Figure 8. Isocontour of the lapse function A'^ for the configuration shown in Fig. ^. 
The plots are cross section is the x = 0, y = and z = planes (note that the x 
coordinate is shifted by 20.5 km with respect to that of Fig. M). The dot-dashed line 
denotes the boundary between the inner numerical grid andthe outer compactified 
one (which extends to spatial infinity) , for the grid system centered on the intersection 
between the rotation axis and the orbital plane. 



5. Innermost stable circular orbit 

An important parameter for the detection of a gravitational wave signal from 
coalescing binaries is the location of the innermost stable circular orbit (ISCO), if 
any. In Table ^, we recall what is known about the existence of an ISCO for extended 
fluid bodies. The case of two point masses is discussed in details in Ref. 1^. 

For Newtonian binaries, it has been shown that the ISCO is located at a 
minimum of the total energy, as well as total angular momentum, along a sequence 
at constant baryon number and constant circulation (irrotational sequences are such 
sequences). The instability found in this way is dynamical. For corotating sequences, 



it is secular instead |33, This turning point method also holds for locating ISCO 
in relativistic corotating binaries |Q. For relativistic irrotational configurations, 
no rigorous theorem has been proven yet about the localization of the ISCO by a 
turning point method. All what can be said is that no turning point is present in 
the irrotational sequences considered in Sect. ^ Fig. ^ shows the variation as the 
orbit shrinks of the ADM mass of the spatial hypersurface t = const (which is a 
measure of the total energy, or equivalently of the binding energy, of the system) 
for the Mb = 1.625 sequence. Clearly, the ADM mass decreases monoticaly. 
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Figure 9. Same as Fig. |^ but for the conformal factor of the spatial metric. 



shift N°i dans le plan orbital 




X [km] 



Figure 10. Shift vector of nonrotating coordinates in the orbital plane, for the 
configuration shown in Figs. m9l 
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Composante xy 




-50 50 
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Figure 11. Same as Fig. g but for the component Kxy of the extrinsic curvature 
tensor. The solid (resp. dashed) lines corresponds to positive (resp. negative) values 

of Kxy 



Model 



Existence of an ISCO References 



Newtonian corotating 
Newtonian irrotational 
GR corotating 
GR irrotational 



ISCO ^ 7 > 2 
ISCO ^ 7 > 2.4 
ISCO 4^ 7 > 5/3 
ISCO exists for 7 = cxD 



0, H 



Table 2. Known results about the existence of an ISCO for extended fluid bodies, 
in terms of the adiabatic index 7. 



without showing any turning point. Figure ^ shows the evolution of the total angular 
momentum along the same sequence. Again there is no turning point. The same 
behaviour holds for the A/b = 1.85 Mq sequence, as shown in Figs. |l^ and |l^. 
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Figure 12. Relative variation of tlie central energy density ec with respect to its 
value at infinite separation e™^ as a function of the coordinate separation d (or of the 
orbital frequency r2/(27r)) for constant baryon mass M-q = 1.625 Mq sequences. The 
solid (resp. dashed) line corresponds to a irrotational (resp. corotating) sequence of 
coalescing neutron star binaries (Figure from Ref. [M). 



Irrotational sequences 

EOS: Y=2 K=0.0332 p„„^"'c' 
0.006 I ' , ' , ' , ' , ' , ' r- 
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Figure 13. Relative variation of the central energy density Sc with respect to its 
value at infinite separation ej."^ as a function of the coordinate separation d Jhr 
constant baryon mass sequences with Mb = 1.625 Mq (solid line, same as in Fig. h3) 
and Mb = 1.85 M© (dashed line). 
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Irrotational sequence Mb=1 .625 M, 

EOS: 71=2 K=0.0332 
1.51 I — ' — , — ' — , — ' — , — ' — , — ' — , — ' — , — ' — , — ' — , — ^ 
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Figure 14. Half of the ADM mass of the binary system as a function of the 
coordinate distance d, along the evolutionnary sequence Mb = 1.625 Mq. 

Irrotational sequence Mb=1 .625 M^^,, 



E0S:y=2 ic=0.0332 

14 I ' , ' , ' , ' , ' , ' , r- 
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Figure 15. Total angular momentum of the binary system as a function of the 
coordinate distance d, along the evolutionary sequence Mb = 1.625 Mq. 



6. Conclusion and perspectives 

We have computed evolutionary sequences of quasiequilibrium irrotational configu- 
rations of binary stars in general relativity. The evolution of the central density of 
each star have been monitored as the orbit shrinks. For a compactification ratio 
M/R = 0.14, the central density remains rather constant (with a slight increase, be- 
low 0.1%) before decreasing. For a higher compactification ratio M/R = 0.17 (i.e. 
stars closer to the maximum mass configuration), a very small density increase (at 
most 0.3%) is observed before the density decrease. It can be thus concluded that no 
substantial compression of the stars is found, which means that no tendency to indi- 
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Irrotational sequence Mg=1.85 M, 

EOS: T=2 K=0.0332 
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Figure 16. Half of the ADM mass of the binary system as a function of the 
coordinate distance d, along the evolutionary sequence Afe = 1-85 Mq. 



Irrotational sequence Mb=1 .85 M, 

EOS: 7=2 K=0.0332 
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Figure 17. Total angular momentum of the binary system as a function of the 
coordinate distance d, along the evolutionary sequence Mb = 1.85 Mq. 



vidually collapse to black hole prior to merger is observed. Moreover, the observed 
density increase remains within the expected error (~ 2%, cf. Sect. 2.1) induced by 
the conformally flat approximation for the 3-metric, so that it cannot be asserted that 
this effect would remain in a complete calculation. 

No turning point has been found in the binding energy or angular momentum 
along evolutionary sequences, which may indicate that these systems do not have any 
innermost stable circular orbit (ISCO). 

All these results have been obtained for a polytropic EOS with the adiabatic index 
7 = 2. We plan to extend them to other values of 7 in the near future. We also plan to 
abandon the conformally flat approximation for the 3-metric and use the full Einstein 
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equations, keeping the helicoidal symmetry in a first stage. 
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